1 Εγκατάσταση βιβλιοθήκης

Εγκατάσταση του πρώτου και κύριου πακέτου (και ενός βοηθητικού)
(Αυτό γίνεται μια φορά σε κάθε Η/Υ - εάν έχει εγκατασταθεί δεν χρειάζεται να δοθεί ξανά η εντολή αυτή)

## CRAN repository:
install.packages("easypackages", repos = "http://cran.us.r-project.org")

2 Φόρτωση βιβλιοθήκης

Φόρτωση του εν λόγω πακέτου

library(easypackages)

3 Εγκατάσταση βιβλιοθηκών

Με την εντολή packages('insert_package_name_here') γίνεται η εγκατάσταση των πακέτων που μας ενδιαφέρουν

Δοκιμάστε το

4 Φόρτωση βιβλιοθηκών

Εντολή για φόρτωση των πακέτων

libraries("betapart", "maptree", "ape", "ggdendro", "dendextend", "plyr", "zoo", 
    "ggplot2", "vegan", "igraph", "rnetcarto")

5 Working Directory

Με την ακόλουθη εντολή, βλέπουμε σε ποιόν φάκελο δουλεύουμε:

getwd()

6 Αλλαγή του working directory

Ενώ με την εντολή αυτή, αλλάζουμε τον φάκελο στον οποίο δουλεύουμε και αποθηκεύουμε δεδομένα
[π.χ., setwd("E:/") εάν θέλουμε να δουλεύουμε στον σκληρό δίσκο Ε]

setwd("E:/Academic/Lessons/EMB/Labs/Labs")

7 Εισαγωγή αρχείου

Τώρα ας εισάγουμε το αρχείο το οποίο περιέχει τα δεδομένα τα οποία μας ενδιαφέρουν να αναλύσουμε

plants <- read.csv("G:/Academic/Lessons/EMB/Labs/Labs/Ryukyus.csv", h = T, sep = ";", 
    row.names = 1)

dim(plants) 

Όπως θα διαπιστώσετε, σήμερα θα ασχοληθούμε με τα Ryukyus Islands, τα οποία είναι περισσότερο γνωστά ως Nansei Islands (στα Ιαπωνικά: Nansei-shoto, δηλαδή “Νοτιοδυτικά Νησιά”), ανήκουν στην Ιαπωνία. Αποτελούν ένα νησιωτικό τόξο, το οποίο βρίσκεται σε γεωγραφικό πλάτος 24^ο - 31.0^οΝ και εκτείνεται σε μήκος 1300 km σε έναν άξονα ΝΔ-ΒΑ, ξεκινώντας από το Kyushu και καταλήγει στην Taiwan. Το νοτιότερο νησί είναι το Yonaguni, ενώ το μεγαλύτερο η Okinawa. Τα μεγαλύτερα σε μέγεθος νησιά είναι κυρίως ορεινά, ενώ τα μικρότερα είναι κοραλλιογενή. Καθώς εντοπίζεται στην ζώνη μετάβασης από το τροπικό στο θερμό εύκρατο κλίμα (δηλαδή, βρίσκεται στην υποτροπική ζώνη), οι κλιματικές συνθήκες είναι ήπιες στο αρχιπέλαγος αυτό, καθ’όλη την διάρκεια του έτους, με μέση θερμοκρασία 15^οC το χειμώνα και 28^οC το καλοκαίρι. Βέβαια, υπάρχουν διαβαθμίσεις: στο Β τμήμα του αρχιπελάγους, το κλίμα χαρακτηρίζεται ως υγρό υποτροπικό, ενώ στο Ν τμήμα του ως τροπικό (ταξινόμηση κατά Koppen). Εξαιτίας της γεωγραφικής τους θέσης, στα νησιά αυτά δεν υπάρχει ξηρή περίοδος και τα περισσότερα εξ αυτών δέχονται περισσότερα από 2000 m βροχής, λόγω των μουσώνων και των καλοκαιρινών τυφώνων. Ως εκ τούτου, στα Ryukyus Islands δεν εμφανίζεται η υποτροπική ξηρή ζώνη, η οποία διαχωρίζει τα εύκρατα δάση από τις υγρές τροπικές περιοχές σε άλλα μέρη του πλανήτη και η φυσική τους βλάστηση αποτελείται από πλατύφυλλα αειθαλή δαση (Good, 1845; Maekawa, 1974; Kira, 1991).

Γεωλογικά, το αρχιπέλαγος διαιρείται σε 3 τμήματα:
1. Northern Ryukyus (τα νησιά βορείως του Tokara Strait)
2. Central Ryukyus (τα νησιά μεταξύ του Tokara Strait και του Kerama Gap)
3. Southern Ryukyus (τα νησιά νοτίως του Kerama Gap)

Εικόνα 1. Γεωτεκτονική διαίρεση του αρχιπελάγους Ryukyus.

Τα Ryukyus είναι ηπειρωτικά νησιά και αποτελούν τις κορυφές της κορδιλλιέρας του αρχιπελάγους Ryukyus και έχουν δημιουργηθεί ως αποτέλεσμα της καταβύθισης της τεκτονικής πλάκας των Φιλιππίνων (Ota, 1998).

8 Κατασκευή μήτρας ανομοιότητας

Για τη διερεύνηση των φυτογεωγραφικών σχέσεων διαφόρων περιοχών η κλασσική μέθοδος ανάλυσης συστάδων (cluster analysis), είναι μια από τις πιο ευρέως διαδεδομένες και περισσότερο χρησιμοποιημένες μεθόδους ανάλυσης, όσον αφορά τη μελέτη βιογεωγραφικών σχέσεων (Kreft & Jetz, 2010; Troia et al., 2012; Linder et al., 2012; Saiz et al., 2013). Σύμφωνα επίσης με τους Georghiou & Delipetrou (2010), η ανάλυση της κατανομής των ενδημικών taxa μέσω συντελεστών ομοιότητας (similarity coefficient) και ιεραρχικής ταξινόμησης (hierarchical classification), μπορεί να αποδειθεί ένα χρήσιμο εργαλείο για τη μελέτη παλαιογεωγραφικών και βιογεωγραφικών φαινομένων εντός του ελληνικού χώρου.

Στο σημείο αυτό θέλουμε να δούμε πόσο μοιάζουν (ή όχι) οι κοινότητες των νησιών μεταξύ τους. Για τον σκοπό αυτό, θα δημιουργήσουμε μια μήτρα ανομοιότητας. Υπάρχουν πολλοί τρόποι (π.χ., μέσω της βιβλιοθήκης vegan ή της βιβλιοθήκης betapart) και πολλοί δείκτες (δείκτης Sorensen ή δείκτης Jaccard), ενώ μπορεί να βασιστούμε στην α-ποικιλότητα (αριθμός ειδών) ή την β-ποικιλότητα (εναλλαγή ειδών.). Σήμερα, αρχικά θα κατασκευάσουμε μια μήτρα ανομοιότητας βάσει της β-ποικιλότητας και του δείκτη Jaccard¹, καθώς πρόσφατα (Kreft & Jetz, 2010) αποδείχθηκε ότι για την διαλεύκανση βιογεωγραφικών σχέσεων, είναι πιο αξιόπιστοι οι δείκτες της β-ποικιλότητας σε σχέση με αυτόυς της α-ποικιλότητας.

Πρέπει όμως να επισημάνουμε ότι η ευρεία χρήση της μεθόδου της ανάλυσης συστάδων ανέδειξε αρκετά από τα μειονεκτήματα της (Williams et al., 1971; Everitt, 1979), όπως την παραγωγή συστάδων, ακόμα και όταν αυτές δεν υπάρχουν (Jain et al., 1986), ή το γεγονός ότι η αντικειμένη φύση της ανάλυσης συστάδων αλλοιώνεται εξαιτίας της υποκειμενικής επιλογής της μεθόδου της ανάλυσης και των συντελεστών ομοιότητας, μιας και οι δύο αυτές επιλογές επηρεάζουν το αποτέλεσμα της ανάλυσης (Jackson et al., 1989).

Για αυτόν τον λόγο, στη συνέχεια, θα αναλύσουμε τα δεδομένα μας μέσω μιας άλλης μεθόδου, αυτήν της δικτυακής βιογεωγραφικής ανάλυσης.

Προς το παρόν, ας κατασκευάσουμε μια μήτρα ανομοιότητας:

############## Create a dissimilarity matrix #################

##-------------------------- Alpha diversity ---------------------------------##
# hf.sor <- vegdist(plants, method = 'bray', binary = TRUE) ## Sorensen
# hf.jac <- vegdist(plants, method = 'jaccard', binary = TRUE) ## Jaccard


## This function (vegdist) always computes a dissimilarity index!

## The higher the value of an island pair, more different the floristic
## assemblages between the two islands


##-------------------------- Beta diversity ----------------------------------##


## First way beta.hf.sor <- betadiver(plants, method = 'sor') ## Sorensen
## beta.hf.jac <- betadiver(plants, method = 'jac') ## Jaccard


## Second way
beta.plants <- beta.pair(plants[, 1:1815], index.family = "jaccard")

## beta.jtu --> Turnover

## beta.jne --> Nestedness

## beta.jac --> Total beta diversity

str(beta.plants)

## List of 3
##  $ beta.jtu: 'dist' num [1:325] 0.5302 0.2283 0.0595 0.0361 0.3452 ...
##   ..- attr(*, "Labels")= chr [1:26] "Takeshima" "Ioujima" "Kuroshima" "Tanegashima" ...
##   ..- attr(*, "Size")= int 26
##   ..- attr(*, "call")= language as.dist.default(m = beta.jtu)
##   ..- attr(*, "Diag")= logi FALSE
##   ..- attr(*, "Upper")= logi FALSE
##  $ beta.jne: 'dist' num [1:325] 0.0107 0.4541 0.7533 0.8082 0.2831 ...
##   ..- attr(*, "Labels")= chr [1:26] "Takeshima" "Ioujima" "Kuroshima" "Tanegashima" ...
##   ..- attr(*, "Size")= int 26
##   ..- attr(*, "call")= language as.dist.default(m = beta.jne)
##   ..- attr(*, "Diag")= logi FALSE
##   ..- attr(*, "Upper")= logi FALSE
##  $ beta.jac: 'dist' num [1:325] 0.541 0.682 0.813 0.844 0.628 ...
##   ..- attr(*, "Labels")= chr [1:26] "Takeshima" "Ioujima" "Kuroshima" "Tanegashima" ...
##   ..- attr(*, "Size")= int 26
##   ..- attr(*, "call")= language as.dist.default(m = beta.jac)
##   ..- attr(*, "Diag")= logi FALSE
##   ..- attr(*, "Upper")= logi FALSE

9 Κατασκευή συστάδας

Με βάση την τρίτη ανομοιότητας (beta.jac) από το προηγούμενο βήμα, θα κατασκευάσουμε μια συστάδα (cluster):

hjac.plants <- hclust(beta.plants$beta.jac, method = "average")

Κάντε το ίδιο και για τις άλλες δύο μήτρες ανομοιότητας και αποθηκεύστε το αποτέλεσμα σε δύο διαφορετικά αρχεία με το όνομα hjne.plants και hjtu.plants για τον εγκιβωτισμό και την αντικατάσταση των ειδών, αντίστοιχα.

10 Βέλτιστος αριθμός ομάδων

Στο σημείο αυτό, μας ενδιαφέρει να δούμε ποιός είναι ο βέλτιστος αριθμός των διαφορετικών ομάδων που υπάρχουν στα δεδομένα μας.

Πώς γίνεται αυτό;

Ο - μέχρι στιγμής - καλύτερος τρόπος είναι να χρησιμοποιήσουμε την εντολή kgs() από την βιβλιοθήκη maptree. Η εντολή αυτή χρησιμοποίει μια μέθοδο η οποία μεγιστοποιεί τις διαφορές μεταξύ των ομάδων, καθώς και τη συνεκτικότητα αυτών, βασιζόμενη στη μήτρα ανομοιότητας που έχουμε κατασκευάσει. Η μικρότερη τιμή της kgs() μας υποδεικνύει τον βέλτιστο αριθμό των ομάδων

kgs_hjac.plants <- kgs(hjac.plants, beta.plants$beta.jac)

kgs_hjac.plants

##         2         3         4         5         6         7         8 
## 27.000000  9.022730  9.603249  8.409222  8.480376  8.800032  9.987913 
##         9        10        11        12        13        14        15 
## 10.570905 11.531190 12.505870 13.481461 14.444968 15.167986 16.148047 
##        16        17        18        19        20        21        22 
## 17.170046 18.082180 19.091765 20.037585 21.039077 22.041738 23.018410 
##        23        24        25 
## 24.015828 25.011172 26.000000

Κάντε το ίδιο και για τις άλλες δύο μήτρες ανομοιότητας και αποθηκεύστε το αποτέλεσμα σε δύο διαφορετικά αρχεία με το όνομα kgs_hjne.plants και kgs_hjtu.plants για τον εγκιβωτισμό και την αντικατάσταση των ειδών, αντίστοιχα.

11 Κατασκευή δενδρογράμματος

Τώρα είμαστε πλέον σε θέση να κατασκευάσουμε το προκύπτον δενδρόγραμμα και να οπτικοποιήσουμε τα αποτελέσματα της ανάλυσης μας:

##------------------------- Preparation  -------------------------------------##
df <- plants ## Create a duplicate of our data

labs <- paste("sta_", 1:26, sep = "") ## New labels
                                      ## 26 is the number of islands
                                      ## You can see that with: dim(plants)
                                       

rownames(df) <- labs                  ## Rows should be named as the data in labs
                                      

cut <- 5                              ## Based on KGS

hc <- hjac.plants                     ## It will contain the desired dissimilarity matrix 
                                      


dendr <- dendro_data(hc, type = "rectangle")    ## Convert it to a dendrogram
                                       
  

clust <- cutree(hc, k = cut)          ## Let's see how the clusters are grouped

clust.df <- data.frame(label = names(clust), ## Create a dataframe containing all the data we need
 cluster = clust)
##----------------------------------------------------------------------------##

################################################################################
                #### Don't change anything here ##### 
################################################################################


##----------------------------------------------------------------------------##
height <- unique(dendr$segments$y)[order(unique(dendr$segments$y), decreasing = TRUE)]
cut.height <- mean(c(height[cut], height[cut-1]))
dendr$segments$line <- ifelse(dendr$segments$y == dendr$segments$yend &
   dendr$segments$y > cut.height, 1, 2)
dendr$segments$line <- ifelse(dendr$segments$yend  > cut.height, 1, dendr$segments$line)
##----------------------------------------------------------------------------##


##----------------------------------------------------------------------------##
dendr$segments$clust <- c(-1, diff(dendr$segments$line))
change <- which(dendr$segments$clust == 1)
for (i in 1:cut) dendr$segments$clust[change[i]] = i + 1
dendr$segments$clust <-  ifelse(dendr$segments$line == 1, 1,
             ifelse(dendr$segments$clust == 0, NA, dendr$segments$clust))
dendr$segments$clust <- na.locf(dendr$segments$clust)
clust.df$label <- factor(clust.df$label, levels = levels(dendr$labels$label))
clust.df <- arrange(clust.df, label)
clust.df$cluster <- factor((clust.df$cluster), levels = unique(clust.df$cluster), labels = (1:cut) + 1)
dendr[["labels"]] <- merge(dendr[["labels"]], clust.df, by = "label")
##----------------------------------------------------------------------------##


##----------------------------------------------------------------------------##
n.rle <- rle(dendr$segments$clust)
N <- cumsum(n.rle$lengths)
N <- N[seq(1, length(N), 2)] + 1
N.df <- dendr$segments[N, ]
N.df$clust <- N.df$clust - 1
##----------------------------------------------------------------------------##


##---------------------------    Dendrogram   --------------------------------##
p <- ggplot() + ylab("Beta diversity (Jaccard)") +  ## Change the label 
                                                  
   geom_segment(data = segment(dendr),
      aes(x=x, y=y, xend=xend, yend=yend, size=factor(line), colour=factor(clust)),
      lineend = "square", show.legend = FALSE) +
   scale_colour_manual(values = c("grey60", rainbow(cut))) +
   scale_size_manual(values = c(.1, 1)) +
   geom_text(data = N.df, aes(x = x, y = y, label = factor(clust),  
   colour = factor(clust + 1)),
      hjust = 1.5, show.legend = FALSE) +
   geom_text(data = label(dendr), aes(x, y, label = label),
       hjust = 0, size = 3, show.legend = FALSE) +
   scale_y_reverse(expand = c(0.2, 0)) +
      coord_flip()

p2 <- p + ggtitle("Ryukyus plants") + ## Change the title
                                    
  
theme(plot.title=element_text(size=rel(1.5), lineheight=0.9, 
face="bold", colour="black"), ## Change the font
                             

axis.title.y=element_blank(), axis.text.y=element_blank(),
        axis.ticks.y=element_blank(), 
        axis.title.x=element_text(face="bold", colour="black", size=11)) 

p2

ggsave("Ryukyus plants.pdf", width=20, height=20, ## Save the image
                                                  
 units="cm", dpi=600)
##------------------------------------------------------------------------------

Τα αποτελέσματα μας συμφωνούν με αυτά των Nakamura et al. (2009); Εάν όχι, γιατί πιστεύετε ότι συμβαίνει αυτό;

Εικόνα 2. Αποτελέσματα της ομαδοποίησης των Ryukyus βάσει του δείκτη ομοιότητας Sorensen.

12 Δικτυακή ανάλυση

Όπως αναφέραμε προηγουμένως, η βιογεωγραφική ανάλυση η οποία βασίζεται στην μέθοδο της ανάλυσης συστάδων, εμφανίζει κάποια εγγενή μειονεκτήματα. Ένας τρόπος να τα υπερκεράσουμε, είναι μέσω της δικτυακής βιογεωγραφικής ανάλυσης, όπως πρώτοι πρότειναν οι Carstensen & Olesen (2009). Η μέθοδος αυτή είναι ένα ισχυρό εργαλείο για να κατανοήσουμε πώς συνδέονται μεταξύ τους κάποια νησιά (Dos Santos et al., 2008), καθώς αποτελεί μια ολιστική προσέγγιση του προβλήματος (i.e., της βιογεωγραφικής σύνδεσης κάποιων νησιών), η οποία δεν εμπεριέχει κάποιο βαθμό υποκειμενικότητας και είναι στατιστικώς αποδεδειγμένα καλύτερη από την ανάλυση συστάδων (Proulx et al., 2005). Η δικτυακή βιογεωγραφική ανάλυση είναι μια από τις ταχύτερα αναπτυσσόμενες μεθόδους ανάλυσης βιογεωγραφικών δεδομένων (Carstensen et al., 2012, 2013; Dalsgaard et al., 2014; Kougioumoutzis et al., 2014, 2017) και βασίζεται - ως επί το πλείστον - στο λογισμικό NETCARTO το οποίο αναπτύχθηκε από τους Guimera & Amaral (2005).

Ας φορτώσουμε λοιπόν ξανά τα δεδομένα μας από το Ιαπωνικό αρχιπέλαγος, υπό την απαραίτητη μορφή:

ryukyus <- read.table("G:/Academic/Lessons/EMB/Labs/Labs/Ryukyus.txt", header = TRUE)

head(ryukyus)

Στη συνέχεια, ας τα διαμορφώσουμε κατάλληλα (matrix) για την βιβλιοθήκη igraph:

mat.ryukyus <- matrix(0, max(ryukyus), max(ryukyus))

for (i in 1:nrow(ryukyus)) mat.ryukyus[ryukyus[i, 1], ryukyus[i, 2]] <- ryukyus[i, 
    2]

g.japan <- graph.adjacency(mat.ryukyus, mode = "undirected", weighted = TRUE, 
    diag = TRUE)

Τώρα είμαστε έτοιμοι να τρέξουμε την ανάλυση μας:

b.japan <- netcarto(igraph::get.adjacency(g.japan, sparse = FALSE), seed = 1234, 
    iterfac = 1, coolingfac = 0.995)

Εάν τρέξουμε την ακόλουθη εντολή:

str(b.japan)

## List of 2
##  $ :'data.frame':    1841 obs. of  5 variables:
##   ..$ name         : Factor w/ 1841 levels "1","10","100",..: 21 22 39 40 43 53 54 60 63 64 ...
##   ..$ module       : int [1:1841] 0 0 0 0 0 0 0 0 0 0 ...
##   ..$ connectivity : num [1:1841] -0.0744 -0.0744 -0.0744 -0.0744 -0.0744 ...
##   ..$ participation: num [1:1841] 0 0 0 0 0 0 0 0 0 0 ...
##   ..$ role         : Ord.factor w/ 7 levels "Ultra peripheral"<..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ : num 0.229

Θα δούμε ότι δημιουργήσουμε μια λίστα που περιέχει δύο αρχεία. Το πρώτο μας ενημερώνει για το πώς ομαδοποιήθηκαν τα νησιά που μελετάμε², ενώ το δεύτερο για το modularity (όσο πλησιάζει στη μονάδα, τόσο καλύτερα).

Στη συνέχεια, θα αποθηκεύσουμε το αρχείο αυτό στο Excel:

write.csv(b.japan, "./Japan.csv")

Τέλος, θα ανοίξουμε το αρχείο που δημιουργήσαμε στο Excel και θα διαπιστώσουμε πως έχουν ομαδοποιηθεί τα νησιά που μελετάμε.

13 Εργασία για το σπίτι

Θα χρειαστεί να στείλετε την εργασία σας σε μορφή PDF³ σε αυτό το e-mail και να απαντήσετε στα ακόλουθα ερωτήματα, αφού έχετε περιγράψει το σετ δεδομενων το οποίο έχετε χρησιμοποιήσει:
1. Πόσα είδη έχει;
2. Πόσα νησιά έχει;
3. Ποιό νησί έχει τα περισσότερα και ποιό τα λιγότερα είδη;
4. Με βάση τα αποτελέσματα σας, κατασκευάστε διαφορετικά αρχεία στο excel τα οποία να περιέχουν την μήτρα ανομοιότητας με την μέθοδο Jaccard και την μέθοδο Sorensen για κάθε ένα από τα στοιχεία της β-ποικιλότητας (συνολική ανομοιότητα, αντικατάσταση των ειδών και εναλλαγή των ειδών), καθώς και για την α-ποικιλότητα.
5. Ποιός είναι ο βέλτιστος αριθμός των ομάδων νησιών στο αρχιπέλαγος του Αιγαίου; Όταν χρησιμοποιήσετε διαφορετική μήτρα ανομοιότητας, θα σας δώσει διαφορετικό αριθμό;
6. Δημιουργήστε την έγχρωμη εικόνα με βάση τον κώδικα που έχετε διαθέσιμο. Πώς ομαδοποιούνται τα νησιά του Αιγαίου; Πώς μπορείτε να αιτιολογήσετε τα αποτελέσματα σας;
7. Τρέξτε την ανάλυση σας με την μέθοδο της ανάλυσης συστάδων και της δικτυακής ανάλυσης μόνο για νησιά του κεντρικού Αιγαίου (35 τον αριθμό). Διαφέρουν τα αποτελέσματα σας μεταξύ των δύο μεθόδων;

Εντολη για να φτιάξετε ένα αρχείο Excel για τις μήτρες ανομοιότητας:

core.jtu <- beta.plants$beta.jtu
core.jne <- beta.plants$beta.jne
core.jac <- beta.plants$beta.jac

cjtu <- as.matrix(core.jtu)
cjne <- as.matrix(core.jne)
cjac <- as.matrix(core.jac)
write.csv
write.csv(cjtu, "./core.pair.jac.spatial.turnover.csv")
write.csv(cjne, "./core.pair.jac.nestedness.csv")
write.csv(cjac, "./core.pair.jac.beta.diversity.csv")